
import galaxychop as gchop # Installed and working 
import numpy as np
import matplotlib.pylab as plt
import matplotlib.gridspec as gridspec
import matplotlib as mpl
import pandas as pd


%matplotlib inline

plt.style.use('dark_background')
plt.style.use('seaborn-bright')

COLOR = 'lightgray'

mpl.rcParams['text.color']      = COLOR  # labels, titles, etc colors
mpl.rcParams['axes.labelcolor'] = COLOR  # Numbers in axis colors
mpl.rcParams['xtick.color']     = COLOR  # ticks colors
mpl.rcParams['ytick.color']     = COLOR  # ticks colors
mpl.rc(      'axes',edgecolor   = COLOR) # axis Colors

mpl.rc('font', size=15) # Font Size

mpl.rcParams['figure.figsize'] = [10, 10] # Size of figures


gal  = gchop.read_hdf5("gal394242.h5")
gal2 = gchop.read_hdf5("gal443801.h5")

print(f'SubhaloID 394242: {gal}\nSubhaloID 443801: {gal2}\n')

SubhaloID 394242: <Galaxy stars=37393, dark_matter=155101, gas=80153, potential=True>
SubhaloID 443801: <Galaxy stars=47789, dark_matter=282694, gas=4459, potential=True>


if not gchop.preproc.is_centered_and_aligned(gal):    
    print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal))
    #gal = gchop.preproc.center_and_align(gal)
    gal = gchop.preproc.center_and_align(gal, r_cut=150)
    print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal))

Centered & Aligned: False
Centered & Aligned: True


if not gchop.preproc.is_centered_and_aligned(gal2):    
    print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal2))
    #gal2 = gchop.preproc.center_and_align(gal2)
    gal2 = gchop.preproc.center_and_align(gal2, r_cut=260)
    print("Centered & Aligned:", gchop.preproc.is_centered_and_aligned(gal2))

Centered & Aligned: False
Centered & Aligned: True

########################################

# Packages 

import numpy as np
import matplotlib.pylab as plt
import illustris_python as il
import matplotlib.gridspec as gridspec

########################################

# Getting the data

basePath = '/home/tnguser/sims.TNG/TNG100-1/output/'
snapNum  = 99  # z=0

fields   = ['SubhaloFlag', 'SubhaloHalfmassRadType',
            'SubhaloMassType', 'SubhaloPos', 'SubhaloVel',
            'SubhaloSFRinRad', 'SubhaloSFR','SubhaloStarMetallicity','SubhaloGrNr']

Subhalo  = il.groupcat.loadSubhalos(basePath=basePath,snapNum=snapNum,fields=fields)

FIELDS=['GroupFirstSub','GroupMassType','Group_M_Crit200','Group_R_Crit200']
HHH=il.groupcat.loadHalos(basePath,99,fields=FIELDS)


GFSub = HHH['GroupFirstSub']
Mtype = HHH['GroupMassType']
M200  = HHH['Group_M_Crit200']
R200  = HHH['Group_R_Crit200']


########################################

# Saving the data

gdms=np.array([0,1,4]) # to keep just columns of gas,dm and stars

Data={}
for i in [394242, 443801]:
    a=Subhalo['SubhaloGrNr'][i]
    Data[str(i)]= dict( id          = i, # The subhalo "ID" is the same as the "index"
                        cpos        = Subhalo['SubhaloPos'][i],
                        cvel        = Subhalo['SubhaloVel'][i],
                        rh_type     = Subhalo['SubhaloHalfmassRadType'][i][gdms],
                        M_type      = Subhalo['SubhaloMassType'][i][gdms],
                        sfr_2rh     = Subhalo['SubhaloSFRinRad'][i],     # in 2 * R_half
                        sfr_tot     = Subhalo['SubhaloSFR'][i],          # in 2 * R_half
                        mh          = Subhalo['SubhaloStarMetallicity'][i],
                        Central     = (i in GFSub),
                        M200_parent = M200[a],
                        R200_parent = R200[a])

np.save('Gal_Data.npy', Data)


Data=np.load('Gal_Data.npy',allow_pickle=True).item()


pd.DataFrame(Data)


x, y, z = gal.stars.x.value , gal.stars.y.value , gal.stars.z.value

r= (x**2+y**2+z**2)**0.5

# print(sum(r>70))

gal.stars.to_dataframe()


PLOT=True

LIM=6

if PLOT: # High Resolution
    x, y, z = gal.stars.x.value , gal.stars.y.value , gal.stars.z.value
    
    mask    = (((x**2+y**2+z**2)**0.5)<LIM)
    x, y, z = x[mask], y[mask] , z[mask] 
    
    H, edges=np.histogramdd((x,y,z),bins=50) #bins=(50,45,55)

    x= (edges[0][1:] + edges[0][:-1])/2
    y= (edges[1][1:] + edges[1][:-1])/2
    z= (edges[2][1:] + edges[2][:-1])/2

    x,y,z = np.meshgrid(x,y,z)

    x,y,z=x.ravel(),y.ravel(),z.ravel()

    #------------------------------------------
    x2, y2, z2 = gal2.stars.x.value , gal2.stars.y.value , gal2.stars.z.value
    
    mask       = (((x2**2+y2**2+z2**2)**0.5)<LIM)
    x2, y2, z2 = x2[mask], y2[mask] , z2[mask] 
    

    H2, edges=np.histogramdd((x2,y2,z2),bins=50) #bins=(50,45,55)

    x2= (edges[0][1:] + edges[0][:-1])/2
    y2= (edges[1][1:] + edges[1][:-1])/2
    z2= (edges[2][1:] + edges[2][:-1])/2

    x2,y2,z2 = np.meshgrid(x2,y2,z2)

    x2,y2,z2=x2.ravel(),y2.ravel(),z2.ravel()


if PLOT: # Scatter density version (plotly)
    import plotly.graph_objects as go
    from plotly.subplots import make_subplots

    #------------------
    # Figure

    fig = make_subplots(rows=1, cols=2,
                    specs=[[{'type': 'volume'},{'type': 'volume'}]],
                    print_grid=False,
                    subplot_titles=( ['SubHalo: 394242','SubHalo: 443801']),
                    horizontal_spacing = 0.075
                    )

    #------------------
    # Scatter plot High resolution

    trace1 = go.Volume(name         = '394242',
                       x            = y,
                       y            = x,
                       z            = z,
                       value        = H.ravel(),
                       isomin       = 1,
                       isomax       = H.max(),
                       opacity      = 0.15, # needs to be small to see through all surfaces
                       opacityscale = 'uniform',#""max",
                       surface_count= 25, # needs to be a large number for good volume rendering
                       colorscale   = 'matter',
                       showscale    = False,
#                        colorbar     = dict(thickness    = 20,
#                                            outlinewidth = 0.7,
#                                            title        = "Count particles per bin",
#                                            titleside    = "right",
#                                            tickmode     = "array",
#                                            tickvals     = [], #[0.001,0.999]
#                                            #labelalias   = {1: str(H.max()), 0:str(H.min()) },
#                                            ticks        = "outside",
#                                            x            = -0.2,)  # position of colorbar
                      )
    
    
    #------------------
    # Scatter plot High resolution

    trace2 = go.Volume(name         = '443801',
                       x            = y2,
                       y            = x2,
                       z            = z2,
                       value        = H2.ravel(),
                       isomin       = 1,
                       isomax       = H2.max(),
                       opacity      = 0.15, # needs to be small to see through all surfaces
                       opacityscale = 'uniform',#"max",
                       surface_count= 25, # needs to be a large number for good volume rendering
                       colorscale   = 'matter',
                       colorbar     = dict(thickness    = 20,
                                           outlinewidth = 0.7,
                                           title        = "Count Stars",
                                           titleside    = "right",
                                           tickmode     = "array",
                                           tickvals     = [], #[0.001,0.999]
                                           #labelalias   = {1: str(H.max()), 0:str(H.min()) },
                                           ticks        = "outside",
                                           x            = 1.2,)  # position of colorbar
                      )
    
    fig.add_trace(trace1, # trace= : figure
                  row=1,  # row=   : Position y
                  col=1)  # col=   : Position x
    
    fig.add_trace(trace2, # trace= : figure
                  row=1,  # row=   : Position y
                  col=2)  # col=   : Position x
    
    #------------------
    # Background color (page and axes)

    # https://www.w3schools.com/colors/colors_rgb.asp

    fig['layout']['paper_bgcolor']="rgb(17,17,17)"
        
    set_bgcolor= dict(showbackground  = True,
                      backgroundcolor = "rgb(101, 101, 101)", #"rgb(212, 219, 249)"
                      gridcolor       = "rgb(10, 10, 10)",
                      zeroline        = False)

    fig.update_scenes(xaxis=set_bgcolor,
                      yaxis=set_bgcolor,
                      zaxis=set_bgcolor)
    
    #------------------
    # Font color and size

    fig.update_layout(font_color = 'white', font_size = 16)

    
    fig.update_layout(height=550, width=1000,title_text="Stars in Halos")
    fig.show()
    print()


g = gal.plot(corner=True)

L=250

g.axes[1, 0].set_xlim(-L, L)
g.axes[2, 0].set_xlim(-L, L)
g.axes[2, 1].set_xlim(-L, L)

g.axes[1, 0].set_ylim(-L, L)
g.axes[2, 0].set_ylim(-L, L)
g.axes[2, 1].set_ylim(-L, L)
g.fig.suptitle('SubHalo: 394242')
plt.show()
print()


g = gal2.plot(corner=True)

L=350

g.axes[1, 0].set_xlim(-L, L)
g.axes[2, 0].set_xlim(-L, L)
g.axes[2, 1].set_xlim(-L, L)

g.axes[1, 0].set_ylim(-L, L)
g.axes[2, 0].set_ylim(-L, L)
g.axes[2, 1].set_ylim(-L, L)
g.fig.suptitle('SubHalo: 443801')
plt.show()
print()


decomp  = gchop.models.JThreshold(reassign=True)


components  = decomp.decompose(gal)


print('SubHalo: 394242')
components.describe()

SubHalo: 394242


gal.plot.sdyn_pairplot(labels=components, sdyn_kws={'reassign': True}, corner=True)

<seaborn.axisgrid.PairGrid at 0x141983640>


g = gal.plot.pairplot(labels=components, palette="Pastel2", corner=True)

g.axes[1, 0].set_xlim(-50, 50)
g.axes[2, 0].set_xlim(-50, 50)
g.axes[2, 1].set_xlim(-50, 50)

g.axes[1, 0].set_ylim(-50, 50)
g.axes[2, 0].set_ylim(-50, 50)
g.axes[2, 1].set_ylim(-50, 50)
g.axes[2, 2].set_ylim(-50, 50)
g.fig.suptitle('SubHalo: 394242')
plt.show()
print()


decomp2     = gchop.models.JThreshold(reassign=True)
components2 = decomp2.decompose(gal2)
print('\nSubHalo: 443801')
components2.describe()

SubHalo: 443801


gal2.plot.sdyn_pairplot(labels=components2, sdyn_kws={'reassign': True}, corner=True)
print()


g = gal2.plot.pairplot(labels=components2, palette="Pastel2", corner=True) 

g.axes[1, 0].set_xlim(-50, 50)
g.axes[2, 0].set_xlim(-50, 50)
g.axes[2, 1].set_xlim(-50, 50)

g.axes[1, 0].set_ylim(-50, 50)
g.axes[2, 0].set_ylim(-50, 50)
g.axes[2, 1].set_ylim(-50, 50)
g.axes[2, 2].set_ylim(-50, 50)
g.fig.suptitle('SubHalo: 443801')
plt.show()
print()

import pandas as pd
df = pd.read_pickle("age_gal394242.pkl")


import pandas as pd
df = pd.read_pickle("age_gal394242.pkl")
age1=df.to_numpy().T[0]

df = pd.read_pickle("age_gal443801.pkl")
age2=df.to_numpy().T[0]

del df


eps1=gal.stellar_dynamics().eps
eps2=gal2.stellar_dynamics().eps


from sklearn.neighbors import KernelDensity # I used a KDE to improve the view of the distirbution

A_plot = np.linspace(-0.2, 14 , 2000)[:, np.newaxis]
e_plot = np.linspace(-1.1, 1.1, 2000)[:, np.newaxis]


A1      = age1[~np.isnan(eps1)][:, np.newaxis]
kde     = KernelDensity(kernel="cosine", bandwidth=0.033).fit(A1)
A1_dens = np.exp(kde.score_samples(A_plot))

A2      = age2[~np.isnan(eps2)][:, np.newaxis]
kde     = KernelDensity(kernel="cosine", bandwidth=0.033).fit(A2)
A2_dens = np.exp(kde.score_samples(A_plot))

e1      = eps1[~np.isnan(eps1)][:, np.newaxis]
kde     = KernelDensity(kernel="cosine", bandwidth=0.033).fit(e1)
e1_dens = np.exp(kde.score_samples(e_plot))

e2      = eps2[~np.isnan(eps2)][:, np.newaxis]
kde     = KernelDensity(kernel="cosine", bandwidth=0.033).fit(e2)
e2_dens = np.exp(kde.score_samples(e_plot))


fig= plt.figure(figsize=(20,10))
gs = gridspec.GridSpec( 2,5,                             # Number of axis y,x
                        height_ratios = [0.3,1],             # relatives ratios of heigh
                        width_ratios  = [0.3,1,0.15,1,0.3], # relatives ratio of with
                        left  = 0.1,      # Space to the edge of left   from the nearest axis
                        right = 0.95,     # Space to the edge of right  from the nearest axis
                        bottom= 0.15,     # Space to the edge of bottom from the nearest axis
                        top   = 0.97,     # Space to the edge of top    from the nearest axis
                        wspace= 0.02,      # Space horizontal between each of the axis
                        hspace= 0.02)      # Space vertical   between each of the axis

C1   = fig.add_subplot(gs[0,0])
C2   = fig.add_subplot(gs[0,4])

ax1  = fig.add_subplot(gs[1,1])
ax3  = fig.add_subplot(gs[0,1],sharex=ax1)
ax5  = fig.add_subplot(gs[1,0],sharey=ax1)

ax2  = fig.add_subplot(gs[1,3],sharey=ax1,sharex=ax1)
ax4  = fig.add_subplot(gs[0,3],sharex=ax1)
ax6  = fig.add_subplot(gs[1,4],sharey=ax1)



#------------------------------------------
# Corners saying density
C1.text(0.5,0.5,'Density',verticalalignment='center',horizontalalignment='center',rotation=45)
plt.setp(C1.spines.values(), visible=False)
C1.minorticks_off()
C1.set_xticks([])
C1.set_yticks([])


C2.text(0.5,0.5,'Density',verticalalignment='center',horizontalalignment='center',rotation=-45)
plt.setp(C2.spines.values(), visible=False)
C2.minorticks_off()
C2.set_xticks([])
C2.set_yticks([])

#------------------------------------------
# Gal 394242

ax1.scatter(age1,eps1,marker='.',alpha=0.1)

ax3.plot(A_plot,A1_dens)
ax3.set_title('Gal 394242')

ax5.plot(e1_dens,e_plot)
ax5.invert_xaxis()

#------------------------------------------
# Gal 433801

ax2.scatter(age2,eps2,marker='.',alpha=0.07)

ax4.plot(A_plot,A2_dens)
ax4.set_title('Gal 433801')

ax6.plot(e2_dens,e_plot)

#------------------------------------------
# Erasing extra label in axis

plt.setp(ax1.get_yticklabels(), visible=False)
plt.setp(ax3.get_xticklabels(), visible=False)

plt.setp(ax4.get_xticklabels(), visible=False)
plt.setp(ax6.get_yticklabels(), visible=False)

#------------------------------------------
# Axes labels
fig.supylabel('$J_z$ / $J_{circ}$',fontsize=22)
fig.supxlabel('Age (Gyr)',fontsize=22)
plt.show()
print()

#Initialize the decomposer KMeans:
decomposer = gchop.models.KMeans(n_components=2, reassign=True) #OR Initialize the decomposer GaussianMixture:
decomposer = gchop.models.GaussianMixture(n_components=2, reassign=True)
components = decomposer.decompose(gal).to_dataframe()


decomposer1 = gchop.models.GaussianMixture(n_components=2, reassign=True)
components1 = decomposer1.decompose(gal)

decomposer2 = gchop.models.GaussianMixture(n_components=2, reassign=True)
components22= decomposer2.decompose(gal2)


gal.plot.sdyn_pairplot(labels=components1, sdyn_kws={'reassign': True}, corner=True)
print()


g = gal.plot.pairplot(labels=components1, palette="Pastel2", corner=True)

g.axes[1, 0].set_xlim(-50, 50)
g.axes[2, 0].set_xlim(-50, 50)
g.axes[2, 1].set_xlim(-50, 50)

g.axes[1, 0].set_ylim(-50, 50)
g.axes[2, 0].set_ylim(-50, 50)
g.axes[2, 1].set_ylim(-50, 50)
g.axes[2, 2].set_ylim(-50, 50)
g.fig.suptitle('SubHalo: 394242')
plt.show()
print()


gal2.plot.sdyn_pairplot(labels=components22, sdyn_kws={'reassign': True}, corner=True)
print()


g = gal2.plot.pairplot(labels=components22, palette="Pastel2", corner=True) 

g.axes[1, 0].set_xlim(-50, 50)
g.axes[2, 0].set_xlim(-50, 50)
g.axes[2, 1].set_xlim(-50, 50)

g.axes[1, 0].set_ylim(-50, 50)
g.axes[2, 0].set_ylim(-50, 50)
g.axes[2, 1].set_ylim(-50, 50)
g.axes[2, 2].set_ylim(-50, 50)
g.fig.suptitle('SubHalo: 443801')
plt.show()
print()


print("\n\n\n\nBy Ian Baeza")



By Ian Baeza

	ptype	ptypev	m	x	y	z	vx	vy	vz	softening	potential	kinetic_energy	total_energy	Jx	Jy	Jz
0	stars	0	5.224283e+05	-0.046012	-0.062069	-0.262818	6.362919	22.511274	-5.774254	0.0	-195699.620206	290.293111	-195409.327095	6.274766	-1.937972	-0.640845
1	stars	0	9.745897e+05	-0.155791	-0.125822	-0.183540	21.987698	8.126047	0.272632	0.0	-196176.962277	274.782915	-195902.179362	1.457155	-3.993156	1.500567
2	stars	0	6.935776e+05	-0.250341	-0.162899	-0.206341	-9.982170	7.308657	15.308136	0.0	-195152.120168	193.699612	-194958.420557	-0.985599	5.891984	-3.455738
3	stars	0	1.070959e+06	-0.223984	-0.253224	-0.378401	-9.688729	-5.382013	-1.178408	0.0	-194695.767625	62.113089	-194633.654536	-1.738158	3.402281	-1.247935
4	stars	0	6.013803e+05	-0.068238	-0.198522	-0.285535	17.203435	7.988667	-16.385347	0.0	-195703.942688	314.128285	-195389.814402	5.533893	-6.030286	2.870127
...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...	...
37388	stars	0	9.867470e+05	134.949071	-161.114980	-11.970927	-59.939303	57.887495	120.366415	0.0	-26553.893402	10715.877926	-15838.015476	-18699.865471	-15525.806786	-1845.255888
37389	stars	0	8.290853e+05	15.012682	15.493519	-4.697090	-13.645752	262.439571	314.563208	0.0	-101280.155757	84005.373302	-17274.782455	6106.393180	-4658.342086	4151.342534
37390	stars	0	1.004103e+06	17.379644	56.147464	-2.132304	-41.672564	304.593269	11.333777	0.0	-64131.846872	47321.058282	-16810.788589	1285.848336	-108.118416	7633.531240
37391	stars	0	1.526682e+06	269.301518	-82.400695	23.001600	-88.927482	-50.437606	-26.101538	0.0	-19698.409864	5566.669709	-14131.740154	3310.930548	4983.709501	-20910.610228
37392	stars	0	1.356703e+06	13.433154	-8.906366	3.706972	172.602906	380.189578	-150.006406	0.0	-111482.516414	98418.900084	-13063.616330	-73.340046	2654.893305	6644.409980

	Particles		Deterministic mass
	Size	Fraction	Size	Fraction
Spheroid	12639	0.338013	1.210456e+10	0.323421
Disk	24753	0.661987	2.532209e+10	0.676579

	Particles		Deterministic mass
	Size	Fraction	Size	Fraction
Spheroid	39829	0.833434	3.740648e+10	0.830121
Disk	7960	0.166566	7.655026e+09	0.169879

Numerical Exercise 3¶

(Deadline: 27/10)¶

Prof. M. Celeste Artale¶

Activity¶

1. Install galaxychop. Download the files and scripts available on the NumEx3 folder. The files gal394242.h5, and gal443801.h5 correspond to SubhaloID: 394242 and 443801, respectively, from IllustrisTNG100-1 at z = 0 (snapshot 99)¶

2. Using jupyterlab in the IllustrisTNG website, investigate the galaxy properties of the two galaxies: stellar masses, star formation rate, metallicity, gas mass and virial mass of the dark matter halo. Are these galaxies central or satellite galaxies?¶

3. Using the script Practico.ipynb, visualize the distribution of gas, dm and stellar particles in real space. Discuss the results found for both galaxies¶

4. Following Practico.ipynb, use the JHistogram method for both galaxies and compare the results obtained. Discuss the stellar mass budget in the disk and bulge.¶

5. Investigate the distribution of the stellar ages of each component, and make a plot of the circularity parameter vs. stellar age for the particles. You will need to use the additional files age_gal443801.pkl and age_gal394242.pkl For reading the files you can use:¶

6. Extend your analysis of the stellar components by using one of the other alternative methods available in galaxychop. Examples:¶

	394242	443801
id	394242	443801
cpos	[63794.09, 12026.507, 5045.7017]	[15467.963, 63014.26, 73353.84]
cvel	[194.05467, 263.4645, 114.92622]	[157.85007, -136.74078, -32.544792]
rh_type	[60.96873, 65.7992, 5.139353]	[109.72671, 74.55964, 1.619801]
M_type	[8.251568, 78.41508, 2.5353656]	[0.48876184, 142.92282, 3.0524666]
sfr_2rh	0.355641	0.0
sfr_tot	0.93635	0.0
mh	0.01982	0.024343
Central	False	True
M200_parent	190.210754	139.27948
R200_parent	201.506332	181.628571